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Abstract. Motor proteins are key players in intracellular transport processes and 
biological motion. Theoretical modeling of these systems has been achieved by the use 
of step processes on one-dimensional lattices. After a comprehensive introduction to 
the total asymmetric exclusion process and some analytical tools, we will give a review 
on different lines of research attracted to the aspects of this systems. We will focus 
on the generic properties of a coupling between the exclusion process and Langmuir 
bulk kinetics that induce topological changes in the phase diagram and multi-phase 
coexistence. 

1 Introduction 

The identification of motion as a manifestation of biological life dates back to the 
earliest records of science itself. The Greek physician Erasistratos of Ceos studied 
biological motion on the length scale of muscles already in the 3rd century BC. 
He imagined muscles to function in the way of a piston contracting and relaxing 
from pneumatic origin. It was not until the invention of the microscope in the 
17th century by van Leeuwenhoek that this theory could be devalidated with 
Swammerdams observation that muscles contract at constant volume. 

Concerning biological motion on a microscopic scale, scientists favored con- 
cepts of "living forces" for many centuries until this was finally ruled out by the 
observations of the Scottish botanist Robert Brown in 1828 who found all kind 
of matter to undergo erratic motion in suspensions. A satisfactory explanation 
was provided by Einstein in 1905 by the interaction with thermally fluctuating 
molecules in the surroundings. However, the molecular details remained unknown 
in the fog of low microscope resolution. Modern experimental techniques [1] have 
lately revealed the causes of sub-cellular motion and transport. 

Today we know that every use of our muscles is the collective effort of a 
class of proteins called myosin that "walk" on actin filaments. Generally spoken, 
we refer to all proteins that convert the chemical energy of ATP (adenosine- 
triphospate) in a hydrolysis reaction into mechanical work as molecular motors. 
These motors are highly specialized in their tasks and occur in a large variety: 
ribosomes move along mRNA strands while translating the codons into proteins, 
dynein is responsible for cilia motion and axonal transport, and kinesins play a 
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key role in cytoskclctal traffic and spindle formation (for an overview see [2] and 
references therein). 

While the exact details of the molecular structure and function of motor 
proteins [3] remain a topic of ongoing research, on a different level attention 
was drawn to phenomena that arise out of the collective interaction of many 
motors. Early research along this line was motivated by mRNA translation that 
is managed by ribosomes. Ribosomes are bound to the mRNA strand with one 
subunit and step forward codon by codon. The codon information is translated 
into corresponding amino acids that are taken up from the cytoplasm and as- 
sembled into proteins. To increase the protein synthesis many ribosomes can be 
bound to the same mRNA strand simultaneously. This fact might induce collec- 
tive properties as was first realized by MacDonald [4] who set up a theoretical 
model for the translation of highly expressed mRNA. The importance of effects 
caused by the concerted action of many motors can be deduced from a very 
simple example that has yet drastic consequences: the slow down of ribosomes 
due to steric hindrance caused by another ribosome in front - comparable to an 
intracellular traffic jam that might significantly slow down protein synthesis. 

A theoretical approach to collective phenomena in intracellular traffic will try 
to simplify the processes of molecular motion down to a single step rate rather 
than focus on the chemical or mechanical details on the molecular level of motor 
steps. Then it becomes possible to model and analyze the behavior of several 
motors with the tools of many-body and statistical physics. We will start this 
review with a short introduction on this single step model in Sec. 2 before we 
introduce the total asymmetric exclusion process (TASEP) as a theoretical model 
for intracellular transport. Sec. 3 describes the stationary states and density 
distributions and their phase diagram as a function of boundary conditions. After 
a review on several recent extensions in Sec. 4, we will focus on the competition of 
TASEP and bulk dynamics in Sec. 5. Before concluding, Sec. 6 contains further 
recent developments. 

2 Model and Methods 

In the quest for a theoretical model for the motion of molecular motors the first 
and simplest choice may be the use of a Poisson process. The "Poisson stepper" 
is assumed to be an extensionless object advancing stochastically in discrete 
steps along a one-dimensional periodic lattice. The process is uni-directional as 
the position of the stepper can be described by x(t) = a n(t) with the discrete 
step size a and the random variable n(t) being a sequence of growing integers. 
Step events occur stochastically with a rate r constant in both space and time. 
Consequently the average time between two steps is then given by the dwell 
time t = 1/r and the probability to find the "Poisson stepper" at a position n 
after time t by the Poisson distribution [5]. 

After we have defined a model for the translocation of a single motor, we 
proceed with our original task which aims at the understanding of collective 
properties of many motors. Of course, more elaborate models have been estab- 
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lishcd [6,7] that account for several rate limiting steps - examples are the ATP 
supply or the availability of amino acids for ribosomal mRNA translation. How- 
ever, the very basic "Poisson stepper" is chosen for reasons of simplicity and in 
order to prevent unnecessary molecular details from masking collective effects. 
Still, the validity and limitations of this simplification have to be kept in mind. 

Being supplied with the dynamics of a single motor, a stage for the concerted 
action of many can now be set up. This was first done in a pioneering work by 
MacDonald [4] and is now widely know as the total asymmetric exclusion process 
(TASEP). It consists of an one-dimensional lattice (Fig. 1) with N sites labeled 
by i = 1 , ■ • • , N and with a spacing of a = L /N, where L is the total length 
of the lattice. For convenience, L is often set to 1 and the lattice spacing then 
referred to as e = 1/N. 



particle '"Vacancy 




Fig. 1. Schematic model of TASEP (particles are injected with rate a, move exclusively 
to the right being subject to hard-core exclusion, and are removed with rate (3) 

Particles have an extension of the size of the lattice spacing and are subjected 
to hard core exclusion due to steric hindrance. Therefore the occupation number 
rii of site i can only take the values or 1. Particles on the lattice attempt 
jumps to their right neighboring site with a rate r, which will be set to unity 
in the following. Hereby, a reference time scale is set. The effective frequency of 
jumps can be much smaller than r when attempted jumps are rejected due to an 
already occupied target site. The attempted jump rate to the left is zero, since 
we deal with a total asymmetric exclusion process, in contrast to the asymmetric 
exclusion process or the symmetric exclusion process, where the jump rate to 
the left is non-zero or even equal to the jump rate to the right. 

Unless one uses periodic boundary conditions, specific dynamic rules have 
to be defined at the boundaries, which play a crucial role in the solution of the 
process. Among different other conditions (reflective, open with a blockage) the 
most common type are open boundaries, which we will use as well: at the left 
boundary (i = 1) particles attempt to attach with a rate a, while they detach 
at the right boundary (i = N) with rate (3. This is equivalent to two additional 
sites i — and i = N + 1 at the boundaries, which are connected to the system 
by the bulk dynamics described above, and are constantly set to the density a 
and 1 — j3 respectively. 
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In spite of its simplicity, TASEP shows a wide range of interesting proper- 
ties. Since it was propelled into the scope of statistical physicists, it has become 
a paradigm for non-equilibrium physics. In contrast to equilibrium systems it 
lacks detailed balance but evolves into a steady state where a non-vanishing 
current is maintained between boundaries. Upon varying these boundary condi- 
tions, TASEP was found to exhibit phase transitions which - following general 
theorems [8] - are not even allowed for one-dimensional equilibrium systems in 
the absence of long-range interactions. However, the analysis of non-equilibrium 
systems is considerably complicated by the lack of universal concepts like the 
Boltzmann-Gibbs ensemble theory. Feasible methods exist nevertheless and will 
be explained in the next section. 

3 Density and Current in Stationary States 

In analyzing exclusion processes research can focus on a multitude of different 
properties. The probably most obvious to address is the density and current dis- 
tribution in the stationary state. This is motivated by two reasons. On the one 
hand, one intuitively attributes a strong importance to density information with 
respect to the biological background as e.g. the ribosome density is connected to 
the rate of protein synthesis. On the other hand, promising experimental tech- 
niques can measure motor densities and may allow for validation of theoretical 
models. Of course, quite extensive research has also been attracted to a multi- 
tude of different properties like correlation functions [9,10], relaxation properties 
[11] or super-diffusive spreading of fluctuations [12] which will not be the topic 
of this review. We will focus on analytical methods (supported by numerical 
simulations) that are designed to investigate spatial density distributions in the 
stationary state of the system. To this end we will introduce some basic tools that 
have proven useful in the exploration of TASEP properties. These are based on 
mean-field approximations and reproduce many results that can also be derived 
exactly. We are well aware that this approach neglects correlations as included 
in the exact solutions that have been achieved for the TASEP density profile by 
applying either recursion relations [13] or a quantum Hamilton formalism with 
Bethe ansatz [14]. 

3.1 Quantum Mechanics and Statistic Properties 

As an introduction we will outline some general statistical properties. At any 
given moment, the system can be found in a certain configuration (i made up of 
the occupation numbers at each lattice site. The next occurring stochastic event 
(i.e. the jump of one particle to a neighboring site) will therefore change the sys- 
tem to another configuration //. The transition probability p^^i is independent 
of the way the system had reached the initial configuration. Since there is no 
memory of the system's history, but any transition probability solely depends 
on the preceding state, TASEP is a Markov process. In order to describe the 
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system's evolution, we can then use a master equation for the probability to find 
the system in a certain state. 

^ = £ [ay^P^t) - awP„(t)] , (1) 

where the w^-,^' are the transition rates from one configuration /i to another fi' . 

How can we now translate this general property into a description of TASEP? 
To this end we will use a convenient notation, which applies methods from the 
quantum mechanics toolbox in order to formulate the master equation in terms 
of operators. It was introduced as "quantum Hamiltonian formalism" and allows 
for exact solutions [14]. We introduce operators rij(i), which act as occupation 
number operators, measuring the presence (m — 1) or absence (m = 0) of a 
particle at site i. This results in the Heisenberg equation (for an introduction 
see e.g. [14]) 

jni{t) = ni_i(i)[l - hi(t)\ - hi(t)[l - h i+1 (t)\ , (2) 

where the first term on the right hand side constitutes the jump of a particle 
from the left neighboring site to site % (and thus a particle gain) and the second 
term a jump from site i to the adjacent lattice site on the right (a particle loss). 
Note the intrinsic exclusion constraint in both terms that prevents jump events 
if the destination site is occupied, i.e. the expression in brackets equals zero. If 
one expresses these gains and losses in current, it becomes convenient to use the 
current operator 

3i(t) = hi(t)[l - h i+1 {i)\ . (3) 

This allows to rewrite (2) as a discretized form of a continuity equation with the 
discrete divergence Vji(t) = ji(t) — ji— 

^ + Vj t (i)=0. (4) 

Similar equations for the boundaries are readily derived in the same way. Since 
we are interested in the average density on a certain lattice site, we need to 
compute the time (or ensemble) average of the operators. Equation (2) gives an 
equation of motion for the operator and allows to solve for the time evolution of 
(nj(i)). In executing the ensemble average of (2) two-point correlation functions 
like (hi-i(t)(l—hi(t))) appear on the right hand side. These correlation functions 
again are connected to higher order correlations via their equations of motion. 
The resulting infinite series of correlation functions can be solved exactly for 
special cases only. Generally, one is required to use mean-field approaches. 



3.2 Mean-Field Solution 



The rather blurry term "mean-field theory" is based on the concept of using time 
or space averages (e.g. by neglecting temporal or spatial fluctuations) and has 
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found a wide range of applications in statistical physics (see [15]). In this chapter, 
we will explain its implementation for the TASEP and show a possible solution. 
To point out the use of mean-field theory in TASEP, we look again at the average 
of the operator hi(t). We are interested in the stationary state and therefore 
averaging signifies either a time or an ensemble average, since the system is 
ergodic. Then the average returns the density at the considered site i as Qi = {hi). 
Performing the average over (2), leaves us with the difficulty of the infinite series 
of correlation functions mentioned earlier. The mean field approximation consists 
now in neglecting any correlations by setting e.g. (hihj) = (hi)(hj) (see [16]). In 
our case, we obtain for the current 

(ni(t)(l - h l+1 (t))) = (ni(t))(l - (n i+1 (t))) . (5) 

This allows then, to rewrite (2) in the stationary state (d,Qi(t)/dt — 0) as 

= - Qi) - Qi(l - Q i+ l) . (6) 

Obviously, (6) could easily be solved numerically, since it forms a system of 
N difference equations. However, it is possible to reduce the set of equations 
to arrive at an explicit solution by making two assumptions. First, note that 
the stationary state condition dgi(t)/dt = implies a conservation of current 
throughout the bulk, as can be seen from (6). Ergo, we just have to solve one 
equation out of the set, to determine the stationary, homogeneous current (and 
density). For that purpose, we use the continuum approximation, which turns 
the spatial lattice variable quasi-continuous. This is achieved for a large number 
N of lattice sites on a lattice of normalized length L — 1. The distribution of 
sites then approaches a continuum as e = L/N <C 1 and x — i/N is rescaled to 
the interval < x < 1. Thereby an expansion in powers of e is allowed. 

g(x ±e) = g(x) ± ed x Q(x) + \e 2 d 2 x g(x) + 0(e 3 ) (7) 

Using this expansion in (6) and the corresponding equations for the boundaries, 
results in the following first-order differential equation if we neglect all terms 
with higher orders in e: 

(2g-l)d x g = 0. (8) 

The corresponding boundary conditions are g(0) = a and g(l) = 1-/3. Because 
we have a first order differential equation that needs to satisfy two boundary 
conditions, we are evidently concerned with an over-determined boundary value 
problem. There are three solutions to (8): £>buik(aO = 1/2 docs not satisfy either 
boundary condition (except for the special case a = (3 = 1/2), while g(x) = C 
can satisfy either the left or the right boundary condition, resulting in g a (x) = a 
and gp(x) — 1 — /?, respectively. 

3.3 Phase Diagram and Domain Wall Theory 

To obtain a general solution satisfying the boundary conditions, both solutions 
need to be matched. Consider the case a,(i < 1/2. To meet the boundary con- 
dition at both sides, the global density function g(x) has to be g a (gp) in an 
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environment close to the left (right) boundary. Since both g a and gp are uni- 
form, the two solutions do not intersect. At this point, we have to go beyond 
mean-field theory and assume that at any given time both solutions are valid in 
non-overlapping areas of the system. Where those areas border, they are con- 
nected by a sharp domain wall (DW) at position x w (see Fig. 2) 

/ n _ f Qa for < x < x w , , . 

Q[X) ~\Qp for x w < x < I . W 

From the dynamics of this domain wall [17] we can deduce important informa- 
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Fig. 2. (left) Schematic density distribution for a = .3,/? = .2: in this situation the 
particle current j a exceeds the current jp, thus carrying more particles to the domain 
wall which is then shifted to the left, (right) Phase diagram of TASEP in a, /3 phase 
space shows a low density (LD), high density (HD) and maximal current (MC) phase 



tion about the system. The key point of domain wall theory is the identification 
of particle currents as the cause of domain wall motion. If for example the cur- 
rent j a — Q a (l — Qa) of the left density solution g a is higher than the current in 
the right part of the system (corresponding to a > (3), particles are transported 
faster to the DW from the left end then they can head on to the right. Thus, the 
domain wall is shifted to the left. The system is filled up, until finally the whole 
bulk density has taken the value of gp except for a small boundary layer 1 at 
the system's left boundary. The opposite happens, if jp > j a , which is the case 
for (3 < a. Hence we can empirically state that the boundary with the smaller 
rate acts as a bottleneck and imposes its density distribution on the system. To 
quantify this behavior, one can use a traveling wave solution [17] of the form 
g(x — Vt) to obtain the domain wall velocity as 

V= j0 ~ ja =0-a. (10) 
gti - g a 

To arrive at a phase diagram in a, /3-space we will analyze the bulk density. A 
positive DW velocity is obtained for a < [3 and pushes the DW to the right 

1 The boundary layer is necessary in order to fulfill the both boundary conditions. 
Its extend is finite for small systems, but will vanish in the thermodynamic limit 
N -> oo. 
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side of the system resulting in a bulk density smaller 1/2 called low-density 
phase (LD). On the contrary, (3 < a leads to a high-density phase (HD), as 
illustrated in the phase diagram (Fig. 2). LD and HD phase are connected by a 
first-order transition. Note, that this phase diagram is equally obtained by the 
above mentioned exact methods. 

The transition line a = /3<l/2in the phase diagram requires special treat- 
ment. The velocity of the domain wall yields zero here, but Monte Carlo simu- 
lations show that the DW actually will make random steps to either side. This 
behavior is caused by the stochasticity of the input and output events that are 
described as rate processes. Hence, the DW makes random steps to either side 
with equal probability, being nothing else than the famous random walk in a 
domain with reflecting boundaries. An average over a sufficiently long time will 
therefore result in a homogeneous probability density over space. The stationary 
density profile will just be a linear slope connecting the two boundary conditions: 

q{x) = a + (1 - - a)x . (11) 

Note that mean-field theory neglects fluctuations and hence fails to predict this 
density distribution, but gives the result of a stationary domain wall. However, 
coming back to the picture of the domain wall as a random walker, the lin- 
ear density distribution can readily be made plausible. Interpreting the random 
walk as free diffusion of the domain wall and keeping in mind that the density 
constraints g(0) — a and g(l) = 1 — hold at all times, the linear density pro- 
file is easily derived as solution of the one dimensional diffusion equation with 
boundary conditions. 

Finally, consider an increase of e.g. a to values that exceed 1/2, while > 1/2 
(crossing the boundary from the upper left to the upper right quadrant in the 
phase diagram). In this situation particles are removed sufficiently fast from the 
system and the supply of particles on the left side acts as a bottleneck limiting 
density and current. For a < 1/2 any increase in a results in an increased 
current (think of a highway, where an additional car will result in a higher overall 
traffic). But above a certain value ac = 1/2 any increase in particle input will 
not increase the current further 2 but will cause the current to diminish again 
(as an additional car will further slow down traffic during rush hour) . As a result 
the bulk will keep its current maximum at g — 1/2 and a boundary layer will 
form at the left side of the system to match the boundary condition. This is 
of course nothing else than the bulk solution gbuik with two boundary layers 
and was baptized maximal current phase (MC). It is reached via second-order 
transitions for values a > 1/2 and > 1/2. A more rigorous treatment of this 
behavior can be gained by computing the collective velocity of the particles [17]. 



2 As j = g(l — g) — g — g 2 has a maximum at g — 1/2. This density- current relation 
is a convenient tool to characterize traffic processes and its plot is often referred to 
as fundamental diagram. 
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4 Biologically Motivated Generalizations of TASEP 



The adoption of lattice gas models for biological systems was followed by a vari- 
ety of efforts to fit TASEP to different realistic environments. For that purpose 
some of the simplifying assumptions TASEP is based on had to be questioned. 
For example, experimental observations have shown that ribosomes typically 
cover an area on the mRNA that exceeds the lattice spacing by a multiple. To 
account for this situation the particles in TASEP have to extend over several 
lattice sites. However it was found that this does not change the phase diagram 
quantitatively [18]. Another direction of research went back to the original field 
of MacDonald to elucidate the importance of initiation and prolongation of ri- 
bosomes [19] or formation of mRNA loops to facilitate the back transport of 
ribosomes from the termination site [20]. 

While particle interactions in TASEP is limited to hard-core potentials, even 
a small increase in the interaction radius - as it could be caused by charged 
molecules - leads to qualitative changes in the phase diagram [21]. It was shown 
that lattice gases with short-range repulsive interaction exhibit a density-current 
relation with two local maxima in contrast to simple TASEP that leads to one 
maximum. This behavior results in a qualitative change of the phase diagram 
that is enriched by four more regions one being a minimal-current phase. 

Further work was dedicated to the scenario of the interaction of different 
species of molecular motors that move in opposite directions. This can either 
happen on the same filament when two different particles are able to surpass 
each other with a jump rate that differs from the jump rate of either particle to 
a free site [22] or on two adjacent one-dimensional filaments [23]. In both cases, 
spontaneous symmetry breaking was observed. 

Intracellular transport along cytoskeletal filaments has also served as a source 
of inspiration for driven lattice gas models. While in the TASEP model motors 
can only bind and unbind on the left and the right boundary respectively, cy- 
toskeletal motors are known to detach from the track to the cytoplasm [2] where 
they perform Brownian motion and subsequently reattach to the track. The 
interplay between diffusion in the cytoplasm and directed motion along the fil- 
ament was studied [24] both in open and closed compartments, focussing on 
anomalous drift and diffusion behavior, and on maximal current and traffic jams 
as a function of the motor density. 

In [25] it has been realized that the on-off kinetics may not only give rise to 
quantitative changes in the transport efficiency but also to a novel class of driven 
lattice gas models. It was shown that the interplay between bulk on-off kinetics 
and driven transport results in a stationary phase exhibiting phase separation. 
This was achieved by an appropriate scaling of the on-off rates that ensures that 
particles travel a finite fraction on the lattice even in the limit of large systems. 
Then, particles spend enough time to "feel" the their mutual interaction and, 
eventually, produce collective effects. In the following section, we will review the 
results of these studies [25,26]. 
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The essential features of cytoskeletal transport are the possibility of bulk at- 
tachment and detachment and a finite residence time on the lattice. The latter 
can be understood as a effect of thermal fluctuations that may overcome the 
binding energy of the motors that is only of the order of several k^T . Hence, 
attachment and detachment is a stochastic process whose dynamic rules have to 
be defined. Parmeggiani et al. [25] chose to use Langmuir kinetics (LK) known as 
adsorption-desorption kinetics of particles on a one- or two-dimensional lattice 
coupled to a bulk reservoir [27]. Particles can adsorb at empty sites and desorb 
from occupied sites and microscopic reversibility demands that the kinetic rates 
obey detailed balance leading to an evolution towards an equilibrium steady 
state describable by standard concepts of equilibrium statistical mechanics. In 
this sense the choice of LK is especially tempting as we are now faced with 
the competition of two representatives of both equilibrium and non-equilibrium 
systems. The system - in the following referred to as TASEP/LK - is defined 
as follows: the well-known TASEP is extended with the possibility of particles 
to attach to the filament with rate u>a and to detach from an occupied lattice 
site to the reservoir with rate wd- According to the type of ensemble (canonical, 
grand canonical) the reservoir is either finite or infinite. Here, the reservoir is 
assumed to be infinite and homogeneous throughout space and time. The den- 
sity on a lattice reached in the equilibrium state of LK is only dependent on 
the ratio K — wd/^a and is completely uncorrelated in both space and time 
for neglection of any particle interaction except hard-core. This is justified by 
the assumption that the diffusion in the cytoplasm is fast enough to flatten any 
deviations of the homogeneous reservoir density. The resulting density profile on 
the lattice is homogeneous and given as Langmuir isotherm ql = Kj (K + 1). 




Fig. 3. Schematic model of TASEP/LK: the TASEP is extended by possible particle 
attachment and detachment in the bulk with rate lua,u~d 

If we now consider the combination of TASEP and LK into the model dis- 
played in Fig. 3, attention has to be paid to the different statistical nature of 
both processes. TASEP evolves into a non-equilibrium state carrying a finite 
current. Since particles are conserved in the bulk, the system is very sensitive 
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to the boundary conditions, whereas LK as a equilibrium process is expected to 
be robust to any boundary effects especially for large systems. Combining both 
processes would thus lead to a trivial domination of LK as the bulk rates iva 
and wd that apply to a large number of bulk sites become predominant over the 
rates a and (3 that only act on the two boundary sites. To observe any interesting 
behavior (i.e. real interplay) between the two dynamics, one needs competition. 
A prerequisite for the two processes to compete are comparable jump rates. To 
ensure that the rates are of the same order independently of the system size, an 
appropriate scaling is needed. To this end a ./V-independent global detachment 
rate J?d is introduced, while the local rate per site scales as 

WD = f ■ (12) 

For the attachment one proceeds similarly. What does this scaling signify phys- 
ically? For an explanation it is instructive, to have a look at the time scales 
involved: a particle on the lattice will perform a certain move on an average 
time scale, which is the inverse of that moves rate. Therefore, a particle spends 
an average time t w I/uj-q on the lattice before it detaches. Bearing in mind 
that the TASEP jump rate is set to unity a particles will jump to its adjacent 
site after a typical time of one unit time step. Therefore the particle will travel 
a number Nt = 1/wd of sites before leaving the lattice. Compared to the lattice 
length this corresponds to a fraction of tit = Nt/N = 1/{Nut>). In order to keep 
this fraction finite in the thermodynamic limit, lod needs to scale as defined in 
(12). Only if the fraction nx is finite, a given particle can experience interaction 
with other particles and give rise to collective phenomena [26] . 



5.1 Mean Field Solution of TASEP/LK 

To obtain density and current distributions of the TASEP/LK, we use again a 
mean-field approach and proceed along the lines of Sec. 3.2. 

First of all, we need to account for the Langmuir kinetics. This is done by 
adding the following terms to the Heisenberg equation (2) of the simple TASEP 
to obtain 

^ii{t) = ni-i(t)(l - hi(t)) - hi(t)(l - h i+1 (t)) - u) D hi(t) + lj a (1 - hi(t)) . (13) 

where the first added term captures the detachment events and the later the 
attachment. Neglecting correlations as done before gives for the stationary state 

= g t -i(l - Qi) - Qi{l - Q i+ i) - w Dft + w A (l - Qi) . (14) 

The boundary conditions are not altered by LK. Using again the power series 
expansion (7) and keeping in mind the scaling of the on and off rates u as in 
(12), we obtain the following ODE: 



-dig + d x g{2g - 1) - us B g + w A (l - q) = 



(15) 
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The ratio K = wa/^d between the attachment and detachment rates will prove 
an important parameter in the analysis of this differential equation. Since the 
case K 7^ 1 is considerably complicated in its mathematical analysis wc refer 
the reader to reference [26] and restrain our discussion to the case K = 1 . In the 
thermodynamic limit (e — > 0) and with uoy> = = & the ODE(15) simplifies 
to first order: 

{d X Q-n)(2 e -i) = o. (i6) 

Obviously there are two solutions to this general ODE problem. The homoge- 
neous density ql = 1/2 given by the Langmuir isotherm and the linear slope 
q(x) = fix + C. The value of C is determined by the boundary conditions and 
leads to the two solutions g a {x) = a + fix and Qp{x) = 1 — /3 — fl + fix. The 
complete density profile g(x) is the combination of one or several of the three 
densities above. Depending on how they are matched, we distinguish several 
phases as explained in the following. 

5.2 Phase Diagram and Density Distributions 

The only area in the phase diagram that does not change compared to simple 
TASEP is the upper right quadrant. This is not surprising since the maximal 
current phase is a bulk controlled regime anyway. Therefore the additional bulk 
dynamics with the Langmuir isotherm at q\, = 1 /2 do not result in any changes in 
the density distribution. In this case, non-equilibrium and equilibrium dynamics 
do not compete but cooperate. 

As in TASEP different solutions can be matched in various ways, the simplest 
being the connection by a domain wall between the left solution g a and the 
right solution Qp. Depending on the current distribution, two possibilities have 
to be distinguished. As both solutions are non-homogeneous, the corresponding 
currents j a and jp will be strictly monotonic (Fig. 4 (left)). If the currents 
equal each other inside the system at a position x w , the DW is localized at this 
position, as a displacement to either side would result in a current inequality 
that drives the DW back to x w (see Fig. 4 (left)). Ergo, the TASEP/LK exhibits 
multi-phase existence of low and high density regions (LD-HD phase) in the 
stationary state on all time-scales, opposed to TASEP where this behavior is 
only observed for short observation times. Recently, this DW localization has 
been observed experimentally [28]. 

If the matching of left and right current is not possible inside the system, 
the known LD and HD phases are found. This is the case for one boundary 
condition being considerably larger than the other and is evidently depending 
quantitatively on the slope of the density solutions. This slope is determined by 
the ratio of the TASEP step rate and the bulk interchange rate fl. For large fl 
any density imposed by the boundary relaxes fast against the Langmuir isotherm 
of £>l = 1/2 resulting a steep slope of the density profile. 

This fact allows for the existence of two other phases with multi-regime co- 
existence. We could imagine a scenario in which the boundary imposed density 
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solutions decay fast enough towards the isotherm to enable a three-regime co- 
existence of low density, maximal current and high density (LD-MC-HD phase). 
Furthermore, a combination of a MC phase with a boundary layer on one side 
and a LD or HD region of finite extend at the other boundary can be imagined. 
Not all these phases will be realized for every value of Q. Instead, the phase 
topology of two-dimensional cuts through the a, (3, Sl-phase space changes. An 
example is shown in Fig. 4 (right). 




a a 



Fig. 4. (left) The DW connects the two densities p a and pg (both dashed) and is 
localized at the point where the correspondent currents j a and jp match. Note the 
finite extend of the DW (localization length) that is only produced for Monte-Carlo 
simulations (solid wiggly line) and is not captured by mean-field results (solid line). 
(right) Topological changes in the phase diagram of TASEP/LK for (a) fi = 0.3, (b) 
n = 0.5, (c) {2=1, from [26] 



5.3 Domain Wall Theory 

After we have derived an analytical solution for the density profile based on the 
mean-field differential equation (15), we now have a closer look on the domain 
wall and its stochastic properties. 

As mentioned before, the density profile exhibits a discontinuity at x w that 
is actually a finite size continuous transition between the high and low density 
for systems of finite size. Only upon increasing the system size N — > oo a sharp 
transition between the left and the right solution occurs. However, this is only 
due to the fact that the lattice spacing decreases to e — ► + . So compared to 
the lattice length L, usually normalized to 1, the domain wall is discontinuous, 
while on the length scale of lattice sites it will still have a finite extend. This 
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extend is an intrinsic statistical feature and is usually referred to as localization 
length. 

The domain wall in its random walk behavior can be cither subjected to equal 
rates (unbiased) as in TASEP or the rates of movement to the left/right can be 
different (biased). In general, the rates will not only be different, but also depend 
on the space variable. To begin with, we will show a way how to derive these 
rates by taking into account fluctuations of particle number [17,29]. Consider a 
situation where all events that can change the particle number (a, fi, fi) have 
a typical time scale that is considerably larger than that of jump processes on 
the lattice. In this case, the time between any entry and exit events is so long, 
that the system has enough time to "rearrange" (to reach a temporary steady 
state) in between . Then it is possible to identify jump rates lo\(x) (lo i (x)) for 
DW movement to the left (right) with the overall rate for entry and exit of any 
particle at any site. Specifically, if a particle enters the system, the DW is shifted 
to the left by a distance of sa e/[gp{x w ) — g a (x w )]. Therefore the rate for the 
DW to move one lattice site to the left is an = uj cn t ry /[gj3(x w ) — g a {x w )]. 

If the density distribution g{x) is known analytically, it is possible to calculate 
Gentry, the overall rate of particle entrance, as the sum over all possible entrance 
events 



The first term captures entrance events from the left boundary reservoir and 
the integral accounts for the Langmuir kinetics. The first multiplicand in both 
terms is the attempted rate of a jump, whereas the difference in brackets states 
the probability of the destination site to be vacant. Along the same lines, the 
exit rate is computed. As we know the analytical density distribution as g{x) = 
a + fix + AO(x — x w ) with the Heavyside function and the DW height 3 
A = 1 — a — fi — H, we can execute the integrals to obtain 



Knowing these rates, one can complete the description of the domain wall as a 
random walker by calculating the position dependent jump rates to the left and 
right uo\{x) and uj i (x). The two rates constitute an effective potential displayed 
in Fig. 5. The DW will always be driven to that point £s in the system where 
^Jcxit(xs) — Gentry (^s)- This position is the center of the DW probability density 
in a stochastic picture and the stationary DW position in a mean-field picture. 
It is in agreement with mean field results and yields 



The other quantity of interest which - in contrast to the DW position - cannot 
be determined by mean-field calculations is the localization length. In order to 

3 The term —Ifl accounts for the diminished height that is caused by the Langmuir 
kinetics on the whole length 1 of the system. 




(17) 




(18) 
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Fig. 5. (left) The DW will be localized at the point where the jump rate to the left u>\ 
and right uj r intersect. This can be interpreted as an potential \ Aui\. (right) Variance 
of DW probability distribution over fl for a = ft = 0.01 and a = j3 = 0.1 in a 
system of N = 500, predictions according to [30] (solid) and based on particle number 
fluctuations (dashed) compared to MC simulations (dots) 



compute this quantity, fluctuations have to be taken into account as we will 
show in the following. If we use the notation p(x) for the probability that the 
domain wall is at position x, then the condition for a stationary DW reads in 
the continuum limit 

uj T (x)p(x) = u>i(x + e)p(x + e) . (20) 

Introducing now y(x) = oj\(x)p(x) and approximate y'(x) = \y(x + e) — y(x)\/e 
we obtain: 

y '( x ) + Ny(x)(l-^) = 0. (21) 
The solution then is given by 

P{x) = ^ = -±--exp[-N [ X dx'(l ^% . (22) 



Zui(x) J X0 lui(x) 

where Z accounts for normalization. In general, Z is not available explicitly, 
but it has been shown [30] that the unnormalized probability function can be 
approximated by a Gaussian 

p(x) cx e -C(*-*s) 2 ; (23) 
where C is given by the second order derivative of the exponent in (22) as 



2dx 2 J Xo V u\(x)J 2cji(xs) 



Hence, the variance a = y/l/(2C) of the domain wall can be easily be obtained 
provided that the jump rates u)\(x) and ui I (x) are available. Evans et al. [30] 
have assumed those rates to be u)\^(x) — j a ,p/{Qp — Qo). Kouyos has shown [31] 
that using the rates (17) derived from the fluctuations of particle number, one 
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arrives at more accurate results compared to Monte Carlo simulations (see Fig. 
5). In this case C evaluates to 

2NDA 

a(l - a) + (3(1 - (3) + n ' 1 ' 



As the width of the DW is given by a — 1/V2C, the localization of the DW 
scales with N~ 1/2 . 



6 Conclusions and Outlook 

Much in the same way as MacDonald's pioneering paper [4] on mRNA trans- 
lation, recent work on kincsin motors walking along microtubules has spurred 
progress in noncquilibrium transport phenomena. The ubiquitous exchange of 
material between the cytoplasm and the molecular track, which originally was 
thought to only lead to quantitative modifications of the dynamics [24] , has re- 
cently been identified [25] as the source for qualitatively new phenomena such 
as phase separation. This introduced a completely new class of non-equilibrium 
transport models. 

These lattice gas models are characterized by a scaling of the on-off rates 
with system size which enables competition of driven motion and equilibrium 
Langmuir kinetics. Hereby, a finite residence time on the lattice ensures coopera- 
tive effects to establish multi-phase coexistence, localized shocks and an enriched 
phase behavior compared to prior TASEP results. 

There are now various routes along which one could proceed. The first one 
is to add more realistic features of the molecular motors, such as the fact that 
they are dimcrs [32] or that there is more than one chemo-mechanical state [28]. 
Such investigations are crucial for a quantitative understanding of intracellu- 
lar traffic in various ways. One might ask how robust the features of minimal 
models are with respect to the addition of more molecular details. In the case 
of dimers the answer is far from obvious since the non-equilibrium dynamics of 
dimer adsorption shows rich dynamic behavior with anomalously slow relaxation 
towards the equilibrium state [27] . How this combines with the driven transport 
along the molecular track was recently analyzed thoroughly [32]. While correla- 
tion effects due to the extended nature of dimers invalidate a simple mean-field 
picture it was found that an extended mean-field scheme can be developed which 
quantitatively describes the stationary phases. Surprisingly, the topology of the 
phase diagram and the nature of the phases is similar to the minimal model with 
monomers. The physical origin of this robustness can be traced back to the form 
of the current-density relation which exhibits only a single maximum. 

The second line of research generalizing the minimal model [25] asks for the 
effect of interactions, more than one molecular traffic lane, "road blocks" such 
as microtubule associated proteins and various other kinds of "disorder", bi- 
directional traffic, coupling of driven and diffusive transport and the like on the 
stationary density profiles and the dynamics. In almost every instance it is found 
that this leads to an even richer behavior with new phenomena emerging. 
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For TASEP it is known that isolated defects (slow sites) depending on their 
strength may either give rise to local density perturbations for low particle den- 
sities or yield macroscopic effects for densities close to the carrying capacity [33] . 
The interplay between coupling to the motor reservoir in the cytoplasm and the 
fluctuations in the capacity limit due to disorder along the track gives rise to a 
number of interesting collective effects [34]. 

Coupling two lanes by allowing particle exchange at a constant rate along the 
molecular track also results in novel phenomena. Similar to equilibrium phase 
transitions described by field theories with two coupled order parameters higher 
order critical points may emerge [35]. 

Coupling diffusive and driven transport, the origin of new phenomena is due 
to a competition of different processes of comparable time scales. The qualita- 
tive failure of mean-field theory in some of these systems [36] comes as quite a 
surprise, since mean-field has proven to predict phase diagrams for large systems 
with an astonishing accuracy in the lattice gas models mentioned above. 
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